In [54]:
%matplotlib inline

import datacube
from datacube_stats.statistics import GeoMedian
import sys

sys.path.append('../Scripts')
from dea_datahandling import load_ard
from dea_plotting import display_map
from dea_plotting import rgb
In [55]:
dc = datacube.Datacube(app='Geomedian_composites')
In [56]:
# Set the central latitude and longitude
central_lat = -28.200
central_lon = 153.169

# Set the buffer to load around the central coordinates - max 0.1 deg
buffer = 0.1

# Compute the bounding box for the study area
study_area_lat = (central_lat - buffer, central_lat + buffer)
study_area_lon = (central_lon - buffer, central_lon + buffer)
In [57]:
display_map(x=study_area_lon, y=study_area_lat)
Out[57]:
In [58]:
# Create a query object
query = {
    'x': (central_lon - buffer, central_lon + buffer),
    'y': (central_lat - buffer, central_lat + buffer),
    'time': ('2019-10-01', '2019-12-31'),
    'measurements': ['nbart_green', 'nbart_red', 'nbart_blue'],
    'output_crs': 'EPSG:32756',
    'resolution': (-10, 10),
    'group_by': 'solar_day'
}

# Load available data
ds = load_ard(dc=dc, 
              products=['s2a_ard_granule', 's2b_ard_granule'], 
              **query)

# Print output data
print(ds)
Loading s2a_ard_granule data
    Applying pixel quality mask
Loading s2b_ard_granule data
    Applying pixel quality mask
Combining and sorting data
Masking out invalid values
    Returning 24 observations 
<xarray.Dataset>
Dimensions:      (time: 24, x: 1967, y: 2219)
Coordinates:
  * y            (y) float64 6.892e+06 6.892e+06 6.892e+06 ... 6.87e+06 6.87e+06
  * x            (x) float64 5.068e+05 5.068e+05 ... 5.264e+05 5.264e+05
  * time         (time) datetime64[ns] 2019-10-03T23:52:49.024000 ... 2019-12-02T23:52:39.024000
Data variables:
    nbart_green  (time, y, x) float32 640.0 564.0 547.0 ... 744.0 773.0 763.0
    nbart_red    (time, y, x) float32 831.0 686.0 628.0 ... 825.0 817.0 803.0
    nbart_blue   (time, y, x) float32 554.0 487.0 441.0 ... 576.0 584.0 569.0
Attributes:
    crs:      EPSG:32756
In [59]:
# Set the timesteps to visualise, this
timesteps = [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23]

# Generate RGB plots at each timestep
rgb(ds, index=timesteps)
In [60]:
# Compute geomedian using all observations in the dataset
geomedian = GeoMedian().compute(ds)
In [61]:
print(geomedian)
<xarray.Dataset>
Dimensions:      (x: 1967, y: 2219)
Coordinates:
  * y            (y) float64 6.892e+06 6.892e+06 6.892e+06 ... 6.87e+06 6.87e+06
  * x            (x) float64 5.068e+05 5.068e+05 ... 5.264e+05 5.264e+05
Data variables:
    nbart_green  (y, x) float32 712.48413 648.0389 ... 730.11786 740.0603
    nbart_red    (y, x) float32 876.9375 781.19324 ... 782.2228 784.9043
    nbart_blue   (y, x) float32 610.8502 552.4766 ... 557.187 558.94324
Attributes:
    crs:      EPSG:32756
In [62]:
# Plot the result
rgb(geomedian, size=10)